function element = implicitReturnMapping(i,element,epsilon)
mu = 0.5*element.E/(1+element.nu);
kappa = element.E/(3*(1-2*element.nu));
lambda = element.E*element.nu/...
    (1+element.nu)/(1-2*element.nu);
%compute trial volume strain
evoltrial = epsilon(1)+epsilon(2);
%compute trial deviatoric strain tensor
edtrial(1) = epsilon(1) - evoltrial/3;
edtrial(2) = epsilon(2) - evoltrial/3;
edtrial(3) = epsilon(3); %engineering notation
edtrial(4) = -evoltrial/3;
%compute trial elastic strain
etrial = (edtrial' - element.ep_old(:,i));
%compute trial deviatoric stress
strial(1) = 2*mu*etrial(1);
strial(2) = 2*mu*etrial(2);
strial(3) = mu*etrial(3);
strial(4) = 2*mu*etrial(4);
%compute trial volumetric pressure
ptrial = kappa*evoltrial;
%check plastic admissibility
phi = yieldFunction();

if phi <= 0
    %update sigma
    sigma(1) = strial(1) + ptrial;
    sigma(2) = strial(2) + ptrial;
    sigma(3) = strial(3);
    element.stress(:,i) = sigma;
    %update tangent operator
    C = C_plane_strain(element.E,element.nu);
    element.constitutiveMatrix(:,:,i) = C;
    %internal variables do not change
else
    %return to smooth portion of the cone
    deltaGamma = computeDeltaGamma();
    J2 = strial(1)^2 + strial(2)^2 + ...
        2*strial(3)^2 + strial(4)^2;
    s = (1 - mu*deltaGamma/sqrt(J2))*strial;
    p = ptrial - kappa*element.etabar*deltaGamma;
    %check validity
    if sqrt(J2)-mu*deltaGamma >= 0
        %update elastic strain
        ee(1) = s(1)/(2*mu) + p/(3*kappa);
        ee(2) = s(2)/(2*mu) + p/(3*kappa);
        ee(3) = 2*s(3)/(2*mu);
        ee(4) = s(4)/(2*mu) + p/(3*kappa);
        %update stress
        tre = ee(1) + ee(2) + ee(4);
        sigma(1) = 2*mu*ee(1) + lambda*tre;
        sigma(2) = 2*mu*ee(2) + lambda*tre;
        sigma(3) = mu*ee(3);
        element.stress(:,i) = sigma;
        %update tangent operator
        
        %update internal variable
        
    else
        %return to apex
        
    end
end

function y = yieldFunction()
end
